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Abstract. In medical SPECT imaging, we seek to simultaneously obtain the internal 
radioactive sources and the attenuation map using not only ballistic measurements but 
also first order scattering measurements. The problem is modeled using the radiative 
transfer equation by means of an explicit nonlinear operator that gives the ballistic 
and scattering measurements as a function of the radioactive source and attenuation 
distributions. First, by differentiating this nonlinear operator we obtain a linearized 
inverse problem. Then, under regularity hypothesis for the source distribution and 
attenuation map and considering small attenuations, we rigorously prove that the 
linear operator is invertible and we compute its inverse explicitly. This allows to 
prove local uniqueness for the nonlinear inverse problem. Finally, using the previous 
inversion result for the linear operator, we propose a new type of iterative algorithm 
for simultaneous source and attenuation recovery for SPECT based on Neumann series 
and a Newton-Raphson algorithm. 


1. Introduction 

1.1. Previous results 

Single-Photon Emission Computed Tomography (SPECT) is a nuclear medicine 
tomographic imaging technique based on gamma ray emmision. The idea is to 
deliver into a patient a gamma-emitting radioisotope (typically technetium-99m) that 
is designed to get attached to certain types of cells or tissues, or distribute in certain 
region, which then start to emit gamma rays (for reference see [15], Chapter 2). This 
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radiation can be measured outside of the patient by a rotating gamma camera which 
can identify both the direction and the energy level of the radiation (see Figure 1 left). 
With the information gathered, the goal is to reconstruct the source distribution of the 
radioisotope inside the patient, hence obtaining an image of the desired specihc tissue in 
study or the region of interest. SPECT is widely used in monitoring cancer treatment 
[39] and also in neuropsychiatric imaging studies [29]. 

The mathematical model commonly used to describe the externally measured 
photons is based on the Radiative Transfer Equation (RTE) (see e.g. the survey [3]) and 
it requires at least two physical parameters: the radioactive source distribution / and 
the attenuation map a. The attenuation map represents the capacity of the medium 
to absorb photons and is, given the medical procedure, an unknown function. As we 
explained before, the radioactive source represents the capacity of the medium to radiate 
photons, and is the main function to be obtained from SPECT. 

The Attenuated Radon transform (AtRT) plays a central role in SPECT, and 
particularly in the extensions made in this work. The inversion formula for AtRT with 
known attenuation was obtained independently by Arbuzov, Bukhgeim and Kazantsev 
in 1998 [1] and by Novikov in 2002 [26] deriving an explicit inverse operator. There 
are several generalization for this result, for general geodesics [33], for complex valued 
coefficients [37] or for more general weight functions [6, 7, 8]. There are also invertibility 
and stability results for partial measurements. In [25], injectivity is obtained by 
measuring in an arbitrarily small open set of angles. Stability for the direct and inverse 
problem can be found at [31] and inversion of data in [2]. 

The identification problem stands in the literature for the simultaneous source and 
attenuation reconstruction of the pair (/, a) from the AtRT. This problem has been 
studied for particular cases of attenuation maps: constant attenuation (when the AtRT 
reduces to the exponential Radon transform) in [18, 34] (see also [19] and the references 
therein), radial attenuation in [28] and piecewise constant attenuation in [11]. For 
particular cases of source functions, the problem has been also tackled in several papers 
[1, 4, 6, 23], and the general non-linear case has been studied in [35]. Nevertheless, in the 
general case, examples of non-uniqueness then appear: for the weighted Radon transform 
in [6], for the exponential Radon transform in [34], for the non-linear identihcation 
problem in [35] and for the linearized one in [4]. 

Another approach to study the identihcation problem is by the characterization of 
the AtRT range. In [24] we can End compatibility properties of the range and in [27] 
there is a full characterization of the range. Recently in [4, 35] some local uniqueness 
and stability results are obtained by using linearization and compatibility conditions for 
the range. 

The identification problem has also motivated several numerical studies. In many 
of them [17, 30, 36], the focus is to hrst obtain a good approximation of the attenuation 
map instead of treating (a, /) as a pair, called attenuation correction algorithms. For 
other numerical aspects and reconstructions see for instance [9, 10, 12, 13, 20, 21, 22, 38]. 
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1.2. Our approach: using lower energy scattering data. 

Our main goal is to reconstruct both the attenuation map and radioactive source 
distribution of an unknown object using the SPECT setting. Although this is the 
same objective as in the above mentioned identification problem, we tackle a different 
inverse problem by using additional scattering measurements. 

Indeed, we can assume that some additional information can be gathered by 
measuring scattered photons outside the object in study. Each time a photon scatters, 
it reduces its energy level (see Figure 1), and gamma cameras can discriminate the 
energy level of photons. Therefore, we can measure separately the gamma rays exiting 
the patient that have not scattered (ballistic photons) and the gamma rays exiting the 
patient that have scattered. Particularly, we are interested in measuring photons that 
just have scattered once (hrst order scattering photons). 

Considering scattering effects leads to introduce a new unknown coefficient in 
the model, an scattering coefficient s(x), that will describe the scattering behavior of 
photons inside the object in study. One of our main assumptions will be to suppose 
some relationship between this scattering coefficient and the attenuation coefficient. 

In summary, we can assume that we can gather more information using the same 
standard device and medical procedure used for SPECT and without the addition 
of new technology or other parameters, except for a change in the protocol for the 
measurements. 



Figure 1. Left: classical SPECT consists in the recovery of the source / by assuming 
known attenuation a only using high energy ballistic photons (solid arrows). We 
propose to use available lower energy photon information corresponding mainly to 
scattered photons (dashed arrows) to recover / and a simultaneously. Right: typical 
distribution of photons in SPECT showing different energy levels (ballistic: higher 
energy levels, scattered: lower energy levels). 


There are three main objectives in this work, the first one is to derive an 
inverse problem that consider scattering effects in the standard mathematical model of 
SPECT that describes the behavior of photons in a medium, this by means of suitable 
assumptions that allows us to deal with the gathered information from the ballistic and 
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first order scattering photons. The second goal is to reconstruct both the attenuation 
and source map from the available data. To achieve this goal we study the operator 
that describes the external measurements by means of a linearization process. The third 
objective is to develop an efficient numerical algorithm that use both the ballistic and 
hrst order scattering photon information to reconstruct both the attenuation and source 
map of an unknown object. 

Notice that the new information given by hrst order scattering photons traveling in a 
certain two dimensional plane contains information of the whole three dimensional body, 
because unlike ballistic photons, scattered photons do not necessarily travel straight 
from the source to the gamma camera. Nevertheless, since the extension to the three 
dimensional case is not too different, and in order to simplify notations, we will restrict 
the analysis and numerical simulations of this paper to the two dimensional case. 

The rest of the paper is organized as follows. In Section 2, we introduce the 
main notations and dehnitions and we develop the mathematical model for the inverse 
problem. In Section 3, we state the nonlinear and linearized inverse problems. In 
Section 4, we hrst state the main mathematical results of invertibility of the linearized 
operator (Theorem 4) and local uniqueness of the nonlinear inverse problem (Theorem 
5) and then we present the proofs. Section 5 contains some numerical experiments that 
illustrated the feasibility of the proposed SPECT using lower energy scattered photons 
in the case of some previously known counterexamples for the identihcation problem 
and for other phantoms. 

Natural extensions of this work are: the three dimensional setting, the use of real 
data, the analysis of more general relationships between attenuation and scattering, 
and the exploration of alternative numerical reconstructions techniques. We intend to 
address these points in a forthcoming paper. 

2. Model description 

2.1. Notation and functional framework 

Let us introduce the notation of the sets and functional spaces used in this paper. Let 
5^ = { 6 ' G : | 6 '| = 1 } be the set of directions in R^, and for 9 = ( 6 ^ 1 , ^ 2 ) G 
61,62 G R, let 6 ^ = (— 6 ^ 2 , 61 ) be its 7 r /2 counterclockwise rotation. Let iC be a compact 
set in R^ of non-empty interior and let iP be a compact set in R^ slightly larger than 
K, for simplicity let us consider K = {x E : \x\ < 1},K = {x G R^ : |x| < 2}. For 
n G N, 0 < a < 1, let C"(R”) be the space of real valued a-Holder continuous functions, 
for m G NU {0, cxd} let C™'(R”) be the space of functions with m continuous derivatives, 
let L^(R”) be the space of square integrable functions and let iL®(R"’),s > 0 be the 
classical Sobolev spaces. In these functional spaces we denote by || • ||c“(R"); II • ||c’"(R")) 
II ■ l|L 2 (]a") and || ■ the usual norms, omitting (R”) in the subscript when the 

context is clear. Also let L“(R”) be the space of essentially bounded functions, and 
abusing the notation let || • ||oo denote the norm in L“(R"') and C'°(R"’). For 12 C R"' 
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let L“(f2) and be the corresponding fnnctional subspaces 

consisting of functions with support contained in hi, this may differ from the standard 
notation, but it is a convenient notation for this article. 

For a function / : R"" x —)■ R let 


ll/llc“(R"xSi) 

and 


sup ||/(•,0)||c-(R") 


||/||h»(R"xS1) 




and let C'"(R"' x S^) and if®(R” x 5^) be the smallest Banach spaces with such norms 
that contain the compactly supported smooth functions (these spaces could have more 
succinctly been described as C"(R"' x S^) := C(S'^; C“(R")) and x : = 

but we adopt the notation commonly used in the related literature, 

see e.g. [25]). 


2.2. Integral operators appearing in the model and the inverse problem 

In the modeling and analysis that will be presented in this work there are a number of 
integral operators that play a crucial role. We proceed to provide a generic definition 
of these operators, leaving the discussion of their properties for the next subsection. 

For a function / : R —)■ R we let FT/ denote its classic Hilbert transform (see e.g. 
[14]). If : Rx 5^ —)■ R then Hg denotes the Hilbert transform of g{-, 0) for each 9 G S^. 
An important integral operator in this paper is the weighted Radon transform. 
Definition 1 (Weighted Radon transform). Let / : R^ —)■ R &e a function and 
w : X ^ R be a weight function, the weighted Radon transform of f, with the 

weight w, is defined as, 

Iwf{s,9)= f w{s9'^ + t9,9)f{s9'^ + t9)dt, s G R, 9eS^. 

Jk 

We will consider specific weight functions that will themselves be composed of 
integral operators, like the beam transform. 

Definition 2 (Beam transform). The beam transform of the function a : R^ —)■ R, of 
the point x G R^, in the direction 9 E S^, is defined as 

poo 

{Ba){x,9) = / a{x + t9)dt, a: G R^, 9 E S^. 

Jo 

The weighted Radon Transform with the exponential of the Beam transform as a 
weight is called the attenuated Radon Transform. 

Definition 3 (Attenuated Radon transform (AtRT)). Let a, / : R^ —?• R, then the 
attenuated Radon transform of f, with attenuation a, is defined as 

Raf{s,9)= f /(s0^ + sgR, 9eS\ 

Jr 

When a = 0 this is called the Radon transform of f and it is denoted as Rf{s,9). 
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Remark 1. The attenuated Radon transform is sometimes defined with a different 
parameterization. We denote such alternative definition as Rff{s, 6) and it satisfies 
Rff{s,9) = Raf{—s,9^). This reparameterization does not change the regularity 
properties of the operator, in particular, results proved for Ra also apply to Rf. 

For the attenuated Radon transform there exists the following inversion formula 
(see [26]). 

Theorem 1 (Inverse of the attenuated Radon transform). Let a, / : —)■ R &e 

continuously differentiable with compact support, then the following complex valued 
formula holds pointwise 

fix) = — Re div [ ■ 9,9) d9, 

47r Jsi 

where his, 9) = 4(/ + iH)R^ais,9) (here i = H is the Hilbert transform and I is 

the identity operator). For a and f less regular this inversion formula holds in a weaker 
sense. 

Definition 4. Let a : R^ —)■ R and J : R x —>■ R, we define J~^is, 9) = J(—s, 9) and 

Rf^Jix) = ^ Re div [ ■ 9,9)d9, x e R^ 

dvr J 51 

where his, 9) = 4(J + iH)R-^ais,9). 

In our analysis we will consider a linearization of the attenuated Radon transform. 
Such analysis will require to work with a weight functional of the following form. 

Definition 5. Let : R^ —)■ R, we define the weight w[u,v] : R^ x —)■ R as 

w[u,v]ix,9) = - + dr,a: G R^, 9 e S\ 

J —00 

And the last integral operator that appears in our modeling and analysis is what 
we call the focused transform. 

Definition 6 (Focused transform). For /, a : R^ —)■ R we define M[a, f] : R^ —)■ R, the 
focused transform of the source f with attenuation a, as 

M[a,f]ix)= f f fix + t9)e-^o<^^^+^oWdtd9, x e R^. 

Js^ Jo 

2.3. Properties of the integral operators and elementary estimates 

The previous integral operators can be dehned in different functional spaces with 
different properties. The following result on the continuity of the weighed Radon 
transform (see e.g. [32]) exemplihes the functional setting in which we will consider 
the integral operators. 
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Theorem 2. If 1/2 < a < 1, f e L‘^{K) and w{x, 6) G x 5^) then 

II-^«)/II_H'1/2(Rx51) ^ C'II'W^I|c“(R2x51)||/I|l2(R2)) 

where the constant C depends only on the compact set K. 

In order to deal with the different nature of the functional spaces involved, the 
remainder of this subsection is devoted to recall some classic results and to provide some 
technical lemmas that will clarify the computations done in section 4.2. A relationship 
between Sobolev and Holder spaces is given by the classic Sobolev embedding. 

Theorem 3 (Sobolev embedding). Let s,n be integers, a >0, if {s — a)/n > 1/2 then 

and the inclusion is continuous. Also for s > 1/2 we have the continuous inclusion 

H"(R) C L“(R). 

The products of functions that appear will fall under one of the following two 
lemmas. Their proofs are obtained by direct calculations. 

Lemma 1. Let H C R” and /i, /2 G then /i ■ /2 G and 

ll/l ■ /2||c“(R") < 2||/i||cQ!(]an)||/2||c«(R")- 
Lemma 2. If f & and g G H^fR), s > 1, then f ■ g E H^^^(R) and 

WfaWm/^iK) < C'll5'lkyR)ll/llrfi/2(R)- 

Next, in Lemmas 3 and 4, we recall some of the basic properties of the Radon 
transform and the beam transform. These properties follow from their definitions or 
from results like the projection slice theorem (see e.g. [25]). 

Lemma 3. We have that 

a) If f{x) = 0 for |x| > 1 then Rf{s,9) = 0,V|s| > 1,V6* G S^. This is also true for 
the weighted Radon transform Iwf{s,9). 

b) Iffe C\K) then \Rf{s,9)\ < C'||/||oo,Vs G R,V0 G Rh 

c) If f & H^{K),t > 0, then '79 G the function s h->• Rf{s, 9) G iL*(R) and 

||R/(-,0)||rf‘(R)<C^||/||HW) y0 eS\ 

The constants C above only depend on K. 

Lemma 4. Let a G C^{K) then (recall .A = {x G R^ : |a;| < 2} ), 

a) 9 ■ dxBa{x, 9) = —a(x), Vx G R^, 79 G , 

b) Ba{x, 9) = 0 if X ■ 9 > 2, 

c) Ba{x, 9) = Ra{x ■ 9^, 9) if x ■ 9 < —2 , 
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d) \Ba{x,e)\ < C'||a||oo,Va; e R2,V0 G S\ 

e) IfaEC^{K) then Ba{x^6) E C^iB? X S^) and \\Ba{x^0)\\c<^i^pi2^gi'^ < C\\a\\c^(^^ 2 y 
The constants above only depend on the compact K. 

We include also the following property for the beam transform. 

Lemma 5. Let a G C^{K) fl H^{K), then 

\Ba{x,e)\ < C||a||Hi(R 2 ),Va; G R^V0 G 

Proof. We have 

\Ba{x,9)\ < ||R|a|(-, 0)1100 < C||i?|a|(-, 0)||jyi(R) < C\\ |a| < C"!|a| 

We used Theorem 3, Lemma 3 and the inequality || |/| ||ji^i(r 2 ) < ||/||/i-i(]a 2 ) (e.g. 

| 161 ). □ 

And to conclude, we present some technical lemmas on Holder regularity for 
functionals that will appear as or in weight functions. 

Lemma 6. Let a G C°‘{K) and let k{x,9) = {x,9) G R^ x S^. Then 

k G ^“(R^ X ^1) and 

||^||c“(i?2xsi) < (l + ||a||c“(R2)) , 

where C is a constant depending only in the compact set K. 

Proof. Fix 9 G S^. Since \k{x,9)\ < and 

\k{x,9) -k{y,9)\ = < e\\^<’^^'>\\^\Ba{x,9^) - Ba{x,9^)\, 

we conclude from Lemma 4 that 

^)l|c“(R2) < )ll°o(^x _|_ ||Ha(-,0‘‘')||ca(R2) < + ||a||c°^(R2), 

where the constant C is independent of 9 and only depends on K. □ 

We also have the following weight appearing in the inversion of the attenuated 
Radon Transform. 

Lemma 7. Let f G H^{K), let h{s, 0) = |(/ + iH)R^f {s, 9) and let ip G C°°{[—2, 2]). 
Then V0 G we have s h-)■ G iL^(R) and 

||(^(.)g±'i( (l + ||/||_f/2(R2)) V0 G S^. 

where C depends only on K and the function p. 
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Proof. From Lemma 3, ||i?/(-, 0)11/^2< C| |/| |h 2 (r 2 ), V0 G and since the Hilbert 
transform is a nnitary operator on iL*(R) > 0, 

||h(-, 0 )||h2(r) < C'||/||h2(r 2).V0 G S'h 

We also have from Lemma 3, 

|g±/*M)| = g±R/M)/2 < gCii/iu Vs G R,V0 G S\ 

With the estimates above, and using the Sobolev inequality 11(71100 < C'llfi'llirpR), h is a 
direct calculation to show that the norm of and its second derivative, are 

controlled as prescribed. □ 

The focused transform will play an important role in our model and in the inverse 
problem. The functional framework in which we will work with it is the following. 

Lemma 8. Let a, f E C°‘{K), then M[a, f] G C“(R^) and 

||M[a,/]||c^<Ce^lHI^(l + ||a||c.)||/||c<., 
were C is a constant depending only in the compact set K. 


Proof. Let us recall that 

M[aJ]{x)= f f + 
isi Jo 

Let kt{x,6) = g-Jo fteix) = f{x + t6) and l^(x) = 1 if a; G R and 0 otherwise. 

Then |l/iel |c“(R 2 ) = ||/||c“(r2), Vf > O,V0 G S^, and as in Lemma 6, 

||^t||c“(R2xsi) < (l + ||a||c'^(]a2)) , Vf > 0. 


Hence 


||M[a,/]||oo < sup 

X Js^ Jo 



\\ftd{-)kt{-, 0)||ool^(a: + tO)dtde 


< (l + ||a||c“(R 2 )) ||/||c“(R 2 ) sup 

X JJo 



1 p^{x + t9)dtd9 


< (l + ||a||c“(R2)) WfWc^ 


(R2). 


Similarly, 


\M[aJ]{x) - M[aJ]{y)\ < 



Jo 


< 



Jo 


\fte{x)kt{x,9) - fte{y)kt{y,9)\ 

■ {tp;{x + t9) + 1 ^( 7 / + t9)) dtd9 

X) 

2||/||c“(R2)||fci||c“(ij2xsi)|a: — y\'^ 


■ (l^(a; + t9) + 1 ^( 7 / + t9)) dtd9 

< + ||a||c“)||/||c“k - 7/1", 


where the constant C only depends on the compact set K. 


□ 
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Lemma 9. Let a, f E C°'{K) and let M = M[a,/] G Then w[a, f] and 

w[a, a ■ M] G C"(R x S^) and they satisfy 

/]||c“(r 2 x 51 ) ^ + ||a||c“)||/||c“5 

||w[a, a ■ M]||cq,(^ 2 x 51 ) < + ||a||c“)^||o||c“||/||c“) 

where C is a constant depending only on the compact K. 


Proof. We recall that 


w[a, f]{x,9) 



J- a(^+te+re)dTj^(^^ ^ X G R^ 0 G 


Following the same steps as in the proof of Lemma 8 we can easily obtain that 


/]||c“(R2xsi) < + ||a||c“)||/||c“- 


This also implies 

||ry[a, a ■ M]||c'a(R2x5i) < + ||a||c“)||a ■ Tf||c“, 

which together with Lemma 1 and Lemma 8 concludes the proof. □ 


2 . 4 . Radiative Transfer Equation model and main simplifying hypotheses 

Let s{x,6,6') be a scattering kernel that gives us the distribution according to which 
photons at the spatial point x G R^, coming from direction 6 E are scattered in the 
direction 9' E S^. The equation that we use to model the propagation of photons with 
attenuation a, source / and scattering s is, for all x G R^ and 9 E 

9 ■ Vxu{x,9) + a{x)u{x,9) + f u{x,9)s{x,9,9')d9' = f{x) + f u{x,9')s{x,9',9)d9' 

isi isi 

lim u{x — t9, 9) = 0. 

t —^~ 1“00 

( 1 ) 

The hrst integral term corresponds to the effect of photons that are scattered away from 
the path dehned by (x, 0), the second integral term is the opposite, and represents the 
gamma rays travelling in the spatial point x G R^ coming from any direction that by a 
scattering process take the path dehned by (x, 9). By introducing the total attenuation: 


a'r(x) = a(x) 

then Equation (1) can be rewitten as 



s{x,9,9')d9' 


( 2 ) 


9-V xu{x,9) + aT{x)u{x,9) = f{x) + / u{x,9')s{x,9\9)d9\x E^ ,9 E . (3) 


'SI 
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Let us introduce Ui{x,9) as the intensity of photons that have been scattered i 
times, thus we can decompose the total intensity u as 

OO 

u{x,9) = y^^Ui{x,9), 

i=0 

(for further reference in this decomposition see e.g. [3]), hence Equation (3) becomes 
the system 

9 ■ 'VxUo{x, 9) + aT{x)uo{x, 9) = /(x), Vx G 9 & 

9 ■VxUi{x,9) + aT{x)ui{x,9) = [ s{x,9,9')ui-i{x,9')d9', Vi > 1, x G R^, 6^ G 5^ 

Js^ 

lim Ui{x — t9, 9) = 0, Vi > 0, x G R^, 9 E S^. 

t —^ ~ 1“00 

(4) 

We hrst assume isotropy of the scattering kernel s{x, 9, 9') = s{x, 9-9') i.e. the scattering 
process just depend on the angle at which photons are scattered, and moreover, we 
assume we can separate variables for the scattering kernel 

s{x,9 -9') = s{x)k{9 -9'). (5) 


Secondly, we assume that the function s(x) is proportional to the attenuation map (i.e. 
3C such that Cs{x) = a(x)) and for the angular variable we assume the scattering 
kernel is independent of the scattering angle (i.e. k{9 -9') = 1/2 t: \/9-9' G [0,1]). With 
these assumptions we have for the total attenuation that 


ar(x) = a(x) + / s{x)k{9 ■ 9')d9'= s{x){C + 1) 
Js^ 

thus redehning C = (27r(l + C))“^ the system (4) becomes 
9 ■ VxUo{x, 9) + aT{x)uo{x, 9) = /(x). 


Vx G R^ 0 G 5^ 


9-VxUi{x,9) + aT{x)ui{x,9) = CaT{x) / Ui_i{x,9')d9', Vi > 1,x G R^, 6^ G 5^ 

Js^ 


lim Ui{x — t9, 9) = 0, 

t^+OO 


Vi > 0, X G R^, 9 G S^. 


( 6 ) 


Proposition 1 . If f and a are uniformly line integrable (i.e. 3D > 

D,\/x G R^,^ G S^) then the system (6) has as unique solution 

/ o 

f{x + t9)e-^^ <^T(-+se)ds^^^ 

•OO 

n*(x, 9) = C [ arix + t9) [ n*_i(x + t9, 9')d9'e -V i > 1. 
J —OO J 

Observe that if f,aE C^{K), then they are uniformly line integrable. 


Proof. The solutions Ui{x,9),i > 0 are obtained by direct integration along the 
characteristics (straight lines) in Equation (6). The line integrability condition ensures 
by induction that the resulting ODEs can be solved uniquely. □ 
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3. Inverse problem 


3.1. Measurements and the Inverse Problem 

We assume that the attenuation arix) and the source f{x) are supported in the compact 
set K, representing the patient. For simplicity from now on we omit the subscript in a^, 
i.e. the total attenuation is named a. For all the other quantities we keep the notation 
of the previous sections. 

As measurements we assume that we are able to record UQ{x,(j)), the ballistic 
photons, and Ui{x,(j)), the hrst order scattering photons, as they exit the patient, i.e. 
we assume the knowledge of uq and ui at all points outside the support of a and /. In 
summary, the inverse problem that we will study is the reconstruction of the source and 
attenuation maps f{x) and a{x) from the measurement of the ballistic and first order 
scattering photons exiting the domain K. 

Under the hypotheses leading to Proposition 1, given a source map f{x) and an 
attenuation coefficient a{x), the intensity of ballistic photons uq{x, 6 ) and the intensity 
of first order scattering photons Ui{x, 6) at any point (x, 0) G x is given by 

/•o Q 

6 )= fix + te)e -X e R^ 0 G 
J —oo 

/ O ^ 

a(x + te)M[a, /](x + a{x+se)ds^^^ x G R^ 0 G 

■OO 

where M[a,/](x) = Uo(x, 0')dd'. Therefore, the ballistic and first order scattering 
photons exiting the domain K correspond to respectively, where we define 

Aiix,9)-.= lim Uiix + T9,9), (x, 0) G R x i = 0,1. (7) 

r^+oo 

The inverse problem can be rephrased as the reconstruction of / and a on K from 
knowledge of the Albedo operator 

A[a, /] = (^ 0 , A) = {(A(2:, 9),Ai{x, 9)), (x, 0) G RA A}. 


Let us write the operator A more explicitly. We have 


Aoix,9) 

Ai{x,9) 

M[a,/](x) 



Vx G R^0 G A, 
Vx G R^0 G S\ 
Vx G R^0 G A, 


and we observe that A and Ai are constant along the directed line define by (x, 9), i.e. 
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Ai{x, 9) = Ai{x + T9, 0), Vr G R, z = 0,1. Abusing the notation we write Vs G R, 6^ G 5^, 
^10(5,0) := A(s 6 '-^, 6 '), 

/ OO 

f{te + 

•OO 

= Ra[f]{s,e), 

Vli(s, 9) := Ai{s9^, 9), 

/ OO 

a{t9 + s9^)M[a,f]{t9 + 

•OO 

= C'R4aM[a,/]](s,0). 

Hence we can write the Albedo operator A as 

Vl[a,/] = (i?4/],C'R4aM[a,/]]), (8) 

and the inverse problem we will study is the inversion of the operator 

[a,f]A{Ra[f],CRa[aM[a,f]]). 


3.2. Formal differential of the Albedo operator and the linearized Inverse Problem 

To study of the invertibility of the Albedo operator A near a known source and 
attenuation pair (a,/), supported in K, we formally compute DA[d,f]{-,-) the 
differential of the Albedo operator at (a, /). 

Since A[a, f] = {Ra[f],CRa[aM[a, f]]) the computation of DA reduces to the 
differentiation of the attenuated Radon transform, which is done in [35], and the 
differentiation of M[a, /]. 

Proposition 2. The formal differential of the Albedo operator at (/, a) is 


DA[dJ]{6a,5f) = 




WiN + flsM 

6a + Ri,(6a ■ M) + + flj(a ■ M[o, 6f]) 


( 9 ) 


where 


w[u, n](x, 9) = — 
M = M[d, f] and 
daMSa{x) = - 



^_Bu{x+Te,e)^^^ + T9)dT, xERf9eS\ 


f{x + t9)e-^oA6c+T9)dr / Sa{x + s9)dsdtd9. 


51 ^0 


Proof. In [35] is shown that the formal differential of (a, /) 1 —>■ Ra [/] is Iw[a,f] [(5a] +Ra [df ], 
which readily implies this result. The computation of daMSa{x) is straightforward. □ 
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To study of the operator DA[d, f]{5a, 5f) we consider some preconditioning. Let 
us recall that the reference pair (a, /) and the perturbation {6a, 5f) are all supported in 
iL = {a: G : |x| < 1}, and also recall that Iu]f{s,6) = Raf{s,9) = 0,V|s| > 1 if / is 
supported in K. We will fix (p G C°°([—2, 2]) such that (p(s) = 1 if |s| < 1 (and (p(s) = 0 
if |s| > 2), and let us define x ^ C°^{K) as x(a;) = ^{\x\). The preconditioning of the 
operator DA[d, f]{5a,5f) consist in the following three steps: 

a) Multiplication component-wise by :p(-)- 

b) Left composition component-wise with R^^- 

c) Multiplication component-wise by x(-)- 

As we will see later in Proposition 9, these three steps map continuously x S^) 

into L^(R^), depending only on d. The resulting preconditioned DA[d, f]{6a, 6f) takes 
the following form as a linear operator {L + Q), that we will consider acting on functions 
supported in K. 

Definition 7. Component-wise we multiply DA\fl, f]{6a,6f) by (f, we left-compose with 
RA and then multiply by x to obtain the operator {L Q) defined as 




Q[d,ffi6a,6f) = 


0 


XRa V4[a,a.M] N + (« ' daMdo) + (u • M[d, 6f]) 


This operator, originally defined on functions {6a,6 f) G {K) x {K), will be extended 
to be acting on functions {6a,6f) G L‘^{K) x L^{K) 

The operator {6a,6f) ^ {L -\- Q){6a,6f) is a preconditioned differential of the 
Albedo operator A, therefore it represents the linearization of the originally non-linear 
inverse problem when {6a, 6f) are supported in K. In section 4 we prove the invertibility 
of the operator {L+Q) in the adequate spaces, providing an explicit inverse and therefore 
solving the linearized inverse problem. 


3.3. Frechet differentiability of the modified Albedo operator 

The term in the differential DA[d, f]{6a, 6f) from the previous subsection, 

quickly exemplifies why the calculated differential is only formal and not a Frechet 
differential: in order to have the required regularity of G x S'^) for 

6a G L^(R^), we need w[d, f] G C"(R^ x S^) for a > 1/2, which is not going to be the 
case for a, f E L^(R^). This obstacle can be overcome by considering a modified Albedo 
operator for the measurements, one arising from a model in which the attenuation and 
the source act in a more regularized way. 

Let e > 0 and G R^ : |x| < 1 — e}. Let : L‘^{Kf) -E H‘^{K) be an 

injective continuous linear operator from L‘^{Kfi) into H‘^{K) that satisfy F^^g > 0 in 
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K ii g > t] m K^. For any g G LS{K^) write g^ '■= F^{g) G H‘^{K) and assnme that 
the transport of photons for a sonrce and attenuation map f,a E is instead 

described by the following modihed Radiative Transfer Equation, 

(f)-Vxu{x,4))-\-ae{x)u{x,(f)) = Cae{x) / u{x, (j)')d(j) + fe{x), xgR^j^gF^. (10) 

isi 

In this case, the measurements are represented by the following modihed Albedo 
operator. 

Definition 8. We define the modified Albedo operator as 

AfiaJ] := A[a,,fi] = (R,J/,], C'R^Ja.Mia,,/,]]). (11) 

On one hand, this modihed Albedo operator adds even more assumptions in the 
model of the measurements. On the other hand, it has a better behaved functional 
structure, which can translate into a more robust implementation in applications. The 
functional structure of Ae is the following. 

Proposition 3. The operator Afia, f] is a well defined operator in the spaces 
Ae : L\Ke) X L‘^{Ke) -E X S^) x x S^) 

and is Frechet differentiable at every point {a, f) G L‘^{Ke) x L'^{Ke), with differential 

DAe[a, f]{Sa,6f) = DAfie, fe]{{da)e, {Sf)e) (12) 

where DA[-, •](•, •) is the formal differential from Proposition 2. 

Proof. Since a, f E L‘^{Ke) then from Sobolev embedding Oe, fe G Ff^{K) C C^{K) for 
a > 1 / 2 . Using Theorem 2 and the Lemmas in Section 2 it follows that Rafife] and 
RafitteMfie, fe]] are in x S^). 

For the Frechet diherentiability it is enough to show that the reminder term in the 
hrst order approximation is quadratic. For both components of the Albedo operator Ae 
we have to study the expansion of the following generic term 

Rb+5b[9 + ^g] = Rbig] + Rb[Sg] + in.ib,g][Sb] + wfig] + w^fig] ( 13 ) 

where lUi = Iwi[b,Sb][g] IU 2 = Iw2[b,5b][^g] with weights 

wfib, 6b] = - 1 + B6h) 

W2[h, 5h] = — 1 ). 

For the hrst component of the Albedo operator we have to take (6, Sb, g, 6g) = 
{oe, {6a)e, fe, {6f)e)■ From Theorem 2 , the Lemmas in Section 2 and the dehnition of Fe, 
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we have 

II^i[/<;]IIh1/2(RxS 1) < C'l|wi||c“(RxSi)||/e||L2(R-) 

< + ||ae||c“(ii'))||-B(^a)e||cO(Rx5i)||-B(^a)e||c“(RxSi)||/e||L2(ii:) 

< + \\ae\\c°‘{K))\\{da)e\\cO{K)\\{Sa)e\\c°‘{K)\\fe\\L'^{K) 

< + ||ajH2(R'))||(5a)e||^2(;^)||/e||L2(R') 

< + ||a||L2(^^))||(5a||i2(x,)||/||L2(R-^) 

where the constant C depends on K, ||a||i 2 (R'^) and 11 /||l 2 (;^^). Similarly 

||W^2[(^/)e]||Hl/2(iaxSl) < C\\5a\\L2(^K,)\\^fWi-^iK,)- 

In other words, for the first component of the Albedo operators, the reminder terms in 
the first order approximation are quadratic and therefore it is Frechet differentiable. 

For the second component of the Albedo operator let Mg = M[ag,/g] and notice 

that 

(a + 5ci)g M[(a + 5a)g, (/ + ^/)e] = hg + (5h)g, 
where hg and {5h)e are given by 

he = OgMg + aedaM[ae, fe]{da)e + agM[ag, (5/)e] + {Sa)eMe, 

{6h)e = Og (M[(a + 5a)g,/g] - Mg - cIaM[ag,/g]((5a)g) 

+ Og (M[(a + Sa)e, {df)e] - M[ag, {6f)e]) 

+ {6a)e{M[{a + 5a)g, /g] - Mg) 

+ (5a)gM[(a + 5a)g, ((5/)g] 

and we expect each term in {Sh)e to be quadratic. 

In the generic Equation (13) we have to take {b,6b, g,Sg) = (og, (5a)g, hg, {Sh)e) and 
now we have to control the reminder 

W = RaAWe] + /.[ag,/^g-agMg][(5a).] + + W2[{6h)e]. 

But since 

\\he\\L^{K) < C, 

||hg - aeMe\\L^(K) < C{\\5a\\L2(K,) + \\Sf\\L^{Ke)): 

\\{bh)e\\L2{K) < C'(||<5a||^2(j^^) + \\Sa\\L2{Ke)\\^f\\L2{K,)), 

the same argument as before will show that the x S*^) norm of the reminder 

term W satisfy a quadratic estimate, and therefore the second component of the Albedo 
operator is also Frechet differentiable. □ 
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4. Main Results and Proofs 

In this section we present the main results and proofs about the linearized inverse 
problem, establishing the appropriate framework for the problem and concluding with 
the invertibility of the operator {L + Q) under some assumptions. The main idea is 
straightforward, to prove the invertibility of the linear operator {L + Q) we will show 
that L is invertible and that Q is a relatively small perturbation. In Section 4, we present 
the main results that build up towards the invertibility of {L + Q) and in Section 4.2, 
we present the proof of the main results and the intermediate technical steps. 


4-1. Main Theorems 

First we describe the functional framework in which we study the operator {L + Q). 

Proposition 4. If d G H‘^{K) and f G C°‘{K),a > 1/2, then the operators L and Q 
from Definition 7 are well defined in the following spaces 

L[dJfQ[dJ] : L\K) X L\K) ^ L^{k) x k{k). 


The second step is the invertibility of the operator L, which is the dominating 
component of the operator {L + Q). 

Proposition 5. Let d E H‘^{K),f E C‘^{K),a > 1/2. Let M{x) = M[d, f]{x) and 
assume that |1/M| is bounded in K. Then the operator L[d, f] is left-invertible, with 
left-inverse 


and 



[a, /] : L^{k) X L^k) -E L\k) x L\k), 

L\ = ( ''t* . ^ 

h ) \ S-Xfl^VGla./llVV] )’ 


(14) 

(15) 


l|i-'[“./lll<2 + C(A-,||a||„a)||l/Af||i„l,V)(l + ll/llc"). 

where C{K, ||a||// 2 ) is non-decreasing in the norm of d. 

The condition |1/M| bounded in k is guaranteed to be fulfilled in the following 

case. 

Proposition 6. If d E L°°{k),f E C'^{k), / > 0 and f ^0, then 1/M E L°°{k) and 

( w \ 2/a 

||/||oo,Va;eiF. 

\j\c°‘ J 


where C is a constant that only depends on diam(i^). 

The next step is to show that the operator Q, i.e. the remainder part of the operator 
(L + Q), is relatively small for d small. Observe that this is immediate in the critical 
case a = 0 since then Q = 0. 
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Proposition 7. Let d E f E C°'{K) with a > 1/2 and ||a||// 2 (if{^ 2 ) < D, then 

f]\\c{L'^{K),L‘2{K)) < C{K, D){1 + ||/||cQ,(ia2))||a||iy2(R2). 

Hence, for a small, the operator (L + Q) is a small perturbation of an invertible 
operator, therefore invertible. 

Theorem 4. Let d E f E C°‘{K) with a > 1/2, / > 0 and / ^ 0. Then L~^ 

and Q are well defined linear operators in the Banach spaces 

Q[d,f],L-^[d,f] : L\k) X L\k) -E L\k) x L\k) 

and there exists H > 0 such that the operator (L + Q)[d, f] defined on Lfi{k) x Lfi{k) 
is left invertible for all d E H‘^{K) satisfying ||a||// 2 (]i^ 2 ) < D. The left inverse is given 
by 

OO ^ 

(L + 0)-‘|a,/] = 5^(-(L-'(2)[a,/]) oL-i[a,/] 

k=0 

and 

||(L + 0 )->[a,/]||< 2 ||L-‘||. 

The constant D depends only in the compact K and the norms ||/||oo, ||/||c“- 

For the non-linear inverse problem with the modihed Albedo operator we can 
provide a local identihcation result for some specihc perturbations. Using the notation 
of Subsection 3.3 let X C Lf{Kf) be a set of functions such that ii g E X then 
ll5'e||L2(x) > C*!Ifi'l|L2(iCe) for a constant C > 0 independent oi g E X. The set X 
can be considered closed under scalarization. 

Theorem 5. Let a,ai,f,fi E L‘^{Kf) with / > 0 and / ^ 0. Assume that 
Sa = ai — a E X and df = fi — fEX. Then there exist constants D,k > 0 depending 
only in the compact set K, the operator F^, the set X and the norms ||/e||oo, ||/e||c“ 
such that */||a||L 2 (^^) < D, | |5a| < k, < n and {a, f) produce the same 

modified measurements as (ai,/a), i.e. if Ae[a, f] = Aefii, fi], then a = ai and f = fi. 

Remark 2. The reguirements on the operator and the space X arise from technical 
considerations and it is not easy to characterize all the pairs {F^,X) satisfying the right 
conditions. Nonetheless, we can provide a simple example that satisfies all the conditions 
and that is of interest in applications. 

For the operator Let e > 0 be sufficiently small and let h^ E U“(R^) satisfy 

he > 0, he ^ 0 and he{x) = 0 if |a:| > e. Define Fe{g) = he* g, the convolution between 

he and g, for g E L‘^{Ke). 

For the set X C Lf{Ke). We say that a partition V of Ke is in the family P if 

each P E V is measurable and contains a ball of radius 2e. For a set A let Xa(x) = I if 

X E A and Xa{x) = 0 otherwise. Define 

X = {gE L\Ke) : g{x) = cpXp{x), Cp G R,P G P}. 

p&v 
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This pair {F^,X) satisfy all the conditions required in the definition of the modified 
Albedo operator (Definition 8) and in Theorem 5 above. 

4-2. Proofs 

This section is devoted to prove the previous results and is organized as follows. We 
start by proving Proposition 6 , which is a direct computation. The next two steps 
consist in obtaining estimates for the operators and We conclude with the 

analysis of the operators L and Q. 

Proof of Proposition 6. Since / G C°‘{K) and / > 0, / 7 ^ 0, there exists x ^ K such 
that f{x) = ll/lloo > 0 . Let 

A={xeK, f(x) > ^}. 

We have 


iRx) - f{y)\ <\f\cfix WyeA^ 


<\f\c^ 


^dist(a:, > 


\x — |/|" \/y e A'^ 

1 /a 

= R. 


2|/|c“ 


Hence B{x,R) C A C K and f{x) > G B{x,R). For x G iF 


M[a,/](x)=/ / f{x + te)e-loA^+^^)ds^^^Q 

Js^ Jo 

P poo 

> e-diamWII“ll- / / f(x + t9)dtde 

isi Jo 


^ g-diam{ft)||a||oo 
^ g-diam(h')||a||oo 



2 -/fii Jo 
ttR^ 


'^Bix,R){x + t9)dtd9, 


2 diam(iF) 


Concluding that 

M[dJ]{x) > 


TT 


2 i+ 2 /odiam(iF) 


diam(X)||a| 


2 /a 


I/I 


Vx G K. 


c° 


□ 


We proceed with some intermediate results needed for the estimates on and 


^ w[a,a-]\4]' 
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Proposition 8. Let a, / G C^{K) with a > 1/2, then 

: L^k) ^ X S'), 

and for 6a G L‘^{K), 

ll-^«)[a,/]['^®]llHi/2(Rxsi) < (1 + ||a||c-) ||/||c“||<5a||2,2, 

ll-^«;[a,a-M[a,/]]['^®]llHi/2(RxSi) ^ (1 + ||a||c“)^ 11«! | C“ 11/1 |c“ 11 11L2 • 

Proof. Follows directly from Theorem 2 and Lemma 9. □ 

Proposition 9. Let d G H‘^{K) and J G x S^), then 

llxi?- ^[</9j]||i2(R2) < (l + ||a||^/2(R2)) II J||h1/2(Rx51)- 

Proof. Define J{s,9) = ip{s)J{—s,9^). Let g G with compact support and 

dehne g = X9- Using the expression for in Definition 4 we write 

(,\A-‘|4=^1,9>i.,k,, = (div f Oel®*'!"'""' (e-'‘ffe*j) [x ■ 0,0)d0,g(x) 


1 


= —Re j (e-^He^.L^ (x ■ 9,9)d9, Vg{x) 


L2(R2) 

L2(R2) 


Since x(x) = 0 for |x| > 2 (hence g{x) = 0 for |x| > 2) and (p{x ■ 912) = 1,V6* G if 
|x| < 2, then 

(xR7i[J],^) =;^Re// ^(x-0/2)0e(^“)U,e")(e-"i/e"j)(x-0,0)d0,V^(x)\ 

\ /l2(r 2) 47r \Jsi V / /l2(r2) 

Dehning F{s,9) = ip{s/2)e~^He’^J{s,9) we get 


L2(R2) 

= —Re (I {9- F{s, 9)) 

4vr \Jx.e=s / l2(rxsi) 


1 




L2(RxS1) 


^Re/ [ g{x)0 ■Vgge^^^^^^’^"Ul{x),F{s,9)) 

4 vr \.lx.e=s / l2(rx51) 


concluding that 



< C" 

/a, /" ^(x)e(^“)U.e^)d/(x),F(s,0)\ 

\ / L2(r2) 


\ Jx.e=s / l2(rxsi) 


C 


^(x)0 • 


' x-6=s 


L2(Kx51) 


(16) 
. (17) 
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Let us bound the terms in (16) and (17). For (16), let k{x,9^) = hence 


A^= (d, j ^dl{x),F{s,9) 

^ a:-6l=s ' l2(rxS1) 

= {dJk9{-s,9^),F{s,9))^,^^^^,^ 

— I li^l/2(]ax5l)ll-^lli^l/2(]ax5l)• 

Using Theorem 2 and Lemma 6 this implies 

Ai < (l + 11011/^2(^^2)) ||^||i,2(r2)||F||//i/2(]i^xS'1)- 

To bound (17) we observe that from Lemma 5, 

\9 ■ ■ Va](x,0)| 

< C'e^ll“ll“||a||// 2 (R 2 ), Vx e G 


hence 


^2 


( [ ~g{x)9 ■V.e^^^^^^'^"Ul{x),F{s,9)) 

\Jx-e=s / l 2 (rxs^) 

< Ce‘^ll“ll°°||a||i/2(R2)||i?|^|(s, 6*) 11^2(^x51)117^(3,6')||i2(Rx5i) 

< C'e‘^''“''°°||a||i/2(R2)||^||/^2(R2)||F(s, 6 *)||l 2 (]r_x 51 )- 


In summary, we have obtained the following bound for all g G C'°°(R^) with compact 
support. 


I (x^d^ J ^ 9 ) I < (1 + ||a||//2) ||F’||_H-i/2(]ax5i)llfi'IU2. 

where F{s,9) = (p(s/2)e“^^®’®^i7e^^^’^V('S)'7(—s, 6*-*-). We complete the proof with the 
following estimate that uses Lemma 1 and Lemma 7, 

< (^gCiiaiu ^ ||s||j^2(ia2))^ f ||e'^('’®V(s)'7(— ,0^)||^i/2(R)d0 

< Ce'^»‘«~ (1 + ||a||«.(K.,)“^_ in(-.«-")liy«|R,<i« 

= (1 + ||a||//2(R2)) ||^||^1/2 (RxS1)- 


□ 
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Proposition 10. Let d E H‘^{K) with ||a||j 72 (R 2 ) < D and f E C^{K) with a > 1/2, 
then 


XRa V4[a,/] 


XRa ^P^w[a,a-M[dJ]] 


C{LHk),L^{K)) 


CiL^K),L^K)) 


<C{K,D)\\f\\c^^n^) 

< C(^, Zl)||a||jy2(R2)||/||(^a(R2), 


with C{K,D) a constant only depending on K and D and can be taken non-decreasing 
in D. 

Proof. This is obtained directly from the estimates in Proposition 8 and Proposition 

9. □ 

Proposition 11. Let d E f E C°'{K) with a > 1/2 and ||a||j:/ 2 (j^ 2 ) < D, then 

Q[dJ] : L^{k) X L^{k) -E L^{k) X 

and 

\\Q[^^ f]\\c{L‘2{Kp,L‘2{Kp) — Cik,D)il + \\f\\ck\ra\\H^- 

Proof. Let 5a, 6f G Lf{k), we bonnd the three terms defining the second component of 
Q[d, f]{Sa, 6f). The first term is bonnded directly from Proposition 10, 

WxRd V-^^[a,a-M]['^®]IUhlFt2) < C {k, H) | | O | || | / | | | | 5a | |x,2(r2) . 

The second term (a ■ daMda) satisfies 

\d-daM6a{x)\ = \d{x) [ [ f{x + te)e-^o^^^+^())<i- f Sa{x + se)dsdtde\ 

Js^ Jo Jo 

< gCllalU 

<C'||a||ooe^ll“ 


h||ooli^(2:) / / \f{xt6 )\dt / 15a(a: + s6*)I 

251 Jn Jr 



00 ; I \6 a{xs6 )\lj^{x)lj^{xs6 )dsd6. 

51 Jr 


Since 1 j^{x)lp-{x + s9) = 0 for |s| > diam(iL),a; G R^, compnting the norm in x 
gives 



OOlUllOO / I I 

51 Jr 


\\d-daM6a\\mn^)<CiK,D)\\d 

<C'(iL,Il)||a||oo||/||oo||5a|U2(R2). 

We can bonnd the third term (a ■ M[d, ■]) similarly since 

/* poo 

\d-M[d,6f]{x)\= d{x) / 6 f{x + t9)e-^J°^^(^+-<^)<i^dtde 

Js^ Jo 

< gC'IPIloo11^1 f f \6f{x + t9)\lj^{x + t9)dtd9, 

Js^ Jr 
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hence 

\\d- M[d,Sf]\\L2t^^2) < C{K,D)\\d\U^^ 

These three estimates readily imply the result. □ 


We have all the estimates needed to prove the main results estated in the previous 
subsection. 


Proof of Proposition 4. The fact that Q[d, f] : LS‘{K) x LS‘{K) —)■ L‘^{K) x L‘^{K) is 
established in Proposition 11. The fact that L[a, /] : L‘^{K) x L‘^{K) —)■ LF‘{K) x L‘^{K) 
is a direct consequence of Proposition 10 and the fact that M = M[a, /] G L°^{K) by 
Lemma 8. □ 

Proof of Proposition 7. This is established in Proposition 11. □ 

Proof of Proposition 5. Given g,h ^ L‘^{K), by the dehnition of L~^, 

2 

< Whiml, + iWoWh + 2|lxfls-V4MlV*llli.. 

L2xL2 

By Proposition 10, for D = ||a||j 72 , 



lIxfls'VGMlVMllU. < C(A-,D)||/||o.||ft/A!f|U>. 


hence 

proving the proposition. 


9 

h 


< 2\\g\\l, + C{k,D){l + \\f\\ck"\\l/M\\l^^K) 


L2xL2 


2 

L2) 


□ 


Proof of Theorem 4. From Proposition 5, Proposition 6 and Proposition 7, 

||L-‘Q||<||L-‘||.||(2|| 

< (2 + C(/?,-D)(l + ||/||c.) ) (C(A,-D)(l + Il/Ilc.)l|a||a.) 

\ ll/||oo / 

< D \C{k,D){l + ||/||c.)^(l + 

1 

< - 
“ 2 

for D sufficiently small, since C{K,D) > 0 is non-decreasing in D. Hence (/ + L~^Q) 
is invertible and its inverse can be written as a Neumann series. Therefore {L + Q) is 
left invertible and its left inverse can be written as 

(L + Q)-' = {I + L-^Q)-^L-^ 

oo 

= Y.HL-'Q)fL-'. 

k=0 

From the estimate for ||L“^Q|| we also get ||(L-|-Q)“^|| < 2||L“^||. □ 
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Proof of Theorem 5. Assume that (a,/), (oi,/i) G x satisfy the 

hypothesis of the theorem. In particular assume that Ae{a, f) = Ae{ai,fi). Write 
6 a = ai — a and 5/ = /i — /, from Proposition 3 we have 

A,{a + 6a,f + 6f) = Ae{a, f) + DA[a^, f^]{{6a)^, (5/),) - R 

where R satisfies 


II-R||hi/2(RxS1) < C{\\6a\\L'2{K,) + ||^/||i, 2 (r-^))^. 

Since the value of the Albedo operator at (a, /) and (oi, /i) agree, we have 

DA[a,,m6a),,{6f),) = R. 

In the relationship above we apply the preconditioning steps before Definition 7: 
component-wise multiply by (p(s), compute the inverse of the attenuated Radon 
transform with attenuation a^, then multiply by x{^)- We obtain 

{L + Q)[a,J,]{{6a),,{6f),) = xR"/[(pR], 


hence 


((fa).,W).) = (L + (3)-‘|a.,/.](x-R-‘|4J-R]). 

The hypotheses for Theorem 4 and Proposition 9 are satisfied, therefore 
||(((5a)e, {6f)e)\\L‘2{K)xL2{K) < C'(||<^a||L2(K,) + ||<^/||l2(k,))^- 
Since 6a, 6f G X and | |5a| |x, 2 (r-) < k, ||^/||l 2 (^) < n then 

\\^a\\L-2(K,) + ||^/||l2(r'^) < CK{\\6a\\L-2(K,) + \\6f\\L-2(K,))- 
For K small enough this implies 5a = 0 and 6f = 0. □ 

5. Numerical experiments 

We now present a MatLab implementation of the Newton-Raphson algorithm based on 
the linearized inverse problem. 

The computational domain is the unit square [—1,1]^ discretized into an equispaced 
cartesian grid of size N x N with N = 256. The quantities of interest (a, /) are 
supported inside the unit disc D = {x"^ + < 1}. The computation of the forward 

measurement operator consists in computing the ballistic and single scattering parts 
(resp. v4,o[a,/] and Ai[a, f] as defined in (7), call these measurements A), outgoing 
traces of the solutions uq,ui of system (4). Such a method is referred to as the iterated 
source method (see e.g. [5]) for solving (3), in the exact same way that in system (4), the 
term Mj_i yields a source term for a transport equation satisfied by Uj. Computing such 
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quantities is based on discretizing uniformly and, for each 6 in this discretization, 
integrating hrst-order ODEs along lines of hxed direction 6 . The latter task is done by 
computing rotated versions of the map one desires to integrate (e.g. a or /) so that the 
direction of integration coincides with one of the cartesian axes of the image, and the 
integration along each row is done via cumulated sums, see [5] for details. In the present 
case, computing the values of a rotated image is achieved via bilinear interpolations. 

Subsequently, the iterative inversion is done by implementing the modihed Newton- 
Raphson scheme 


(a“,/“) = (0,l), 

(«”+',/”+■) = (o”,/”) - L-‘ ( ^(-OL-')‘[F,(a”,r)] ) -.4), 


. fc =0 


(18) 


where the operators L, Q come from Dehnition 7 and the modihed albedo operators 
comes from Dehnition 8. 


Remark 3. This algorithm is a modified Newton-Raphson algorithm in the sense that 
we use the inverse of DAlFfia"^, f'^)] instead of the inverse of DAe[Ffia"‘, f"')] in the 
right hand side of eguation (18). The operator F^ is introduced in Section 3 for 
theoretical purposes and it can be chosen to he an approximation of identity (mollifier). 
Numerically, doing so adds robustness to the scheme and does not affect the convergence 
of the algorithm to the correct target functions. 

In all experiments below, 8 iterations of the scheme (18) are enough to ensure 
convergence, and the Neumann series approximated by its hrst 4 

terms. The implementation of L~^ and Q is straightfoward via rotations, cumulated 
sums and pointwise multiplications/division on the cartesian grid. 

Axes on figures. In Figures 2 through 6, functions of {x, y) are represented on the 
unit square [—1,1]^. For f = 0, 1, the measurement data Afia, f] (e.g. on Fig. 2, bottom 
row) are represented by their values for Ai{s9-^, 9) for 6^ G [0, 27r] on the horizontal axis 
and s G [—1,1] on the vertical axis. 

In the sections below, we present two series of experiments. Section 5.1 aims 
at showing that considering the data Ai{a, f) in addition to Ao{a, f) tremendously 
improves the conditioning of an inverse problem referred to as the identihcation problem. 
Secondly, while Section 5.1 treats the reconstruction of smooth pairs (a,/). Section 
5.2 illustrates the performance of our algorithm in the case of discontinuous unknown 
coefficients from measurements with different levels of noise. 


5.1. Non-unigue pairs and trapping geometries 

The problem of reconstructing the pair (a,/) from only the measurements Ao{a, f) 
is known as the identification problem. Recent theoretical work on the identihcation 
problem [35] and corresponding numerical experiments in [20] show that this problem 
is badly conditioned in at least two ways: 
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Lack of injectivity. If a and / are both radial (and smooth enongh), there exists a 
radial fnnction /q snch that Ao{a, f) = .4,o(0,/o). This lack of injectivity prevents 
some experiments done in [20] from converging to the right nnknowns. 

Instability. On the linearized problem (say, reconstrnct {6a, Sf) from <5.4.0 aronnd a 
backgronnd (a,/)), microlocal stability is lost when a certain Hamiltonian flow 
related to the backgronnd (a,/) has trapped integral cnrves inside the domain of 
interest (referred to as a trapping geometry). In this case, experiments done in [20] 
show the presence of artifacts in reconstrnctions. 

The nnmerical experiments of this section aim at showing that accessing the additional 
measnrement Ai{a, f) helps at snccessfnlly reconstrncting both nnknowns (a, /) in both 
scenarios described above. Fignre 2 displays a pair (a, /) corresponding to each scenario, 
as well as the corresponding forward data. 



Figure 2. Top row: a “non-unique” radial pair (ai,/i) (left) and a non-stable pair 
{a 2 ,h) as defined and studied in [20] (( 02 , / 2 ) are rotated by 90 degrees). Bottom row: 
the data (4o(ai,/i),4i(ai,/i)) (left) and (4o(a2,/ 2 ),4i(a2,/ 2 )) (right). 


Pointwise errors on reconstrnctions are shown on Fig. 3, where in both scenarios, 
the reconstrnction is excellent. The cnt plots in the second row of Fig. 3 illnstrate the 
speed of convergence of the method as well as the cross-talk between the reconstrncted 
qnantities dnring the iterations. 


5.2. Reconstruction of non-smooth coefficients and robustness to noise 

We now consider the case of discontinnons nnknown coefficients. We rnn three 
simnlations nsing the same discontinnons nnknowns (shown in Fig. 4), one nsing 
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Figure 3. Top row: pointwise error on (ai,/i) (left) and (021/2) (right) after 
convergence of the Newton-Raphson algorithm. Bottom row: for the reconstruction of 
02 (left) and /2 (right), cut plots of the iterations at {x = 0 }. 

noiseless data and the other two polluted with instrumental noise with different levels. 

Noise model. We add to our measurement a noise of two natures: 

(i) The first kind, modelling instrumental noise, is characterized by an amplitude A so 
that, each data pixel value p is replaced by a draw “y4-Pois(^)”. 

(ii) After this is done, a background noise is added, characterized by a bias value 

/f^added background photons 

H = -. 

#photons measured 

After deciding a value for a quantum q of energy representing one photon, for each 
additional photon, we add g to a pixel chosen at random with uniform probability 
among all data pixels. 

The experiments with “low noise” and “high noise” below are carried out with the 
respective values {A,B) = (0.2, 0.5) and {A,B) = (0.4,5). The forward data 
are displayed on Fig. 5, and the errors after convergence in all three cases (noiseless, 
low noise, high noise) are displayed in Fig. 6. The relative mean square (RMS) errors 
after 8 iterations are summarized in Table 1. 

Comments. 

(i) Regarding the choice of initial guess, / should be chosen at first to be a non¬ 
vanishing function, so as to prevent the vanishing of the focused transform M[a, f] 
which appears in denominators of subsequent operations. As seen above, the choice 
/ = 1 leads to satisfactory convergence. 
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Figure 5. Forward data Ao{a,f) (top row) and Ai{a, f) (bottom row), with (a,/) 
given in Fig. 4. Left to right: noiseless, low noise, high noise. 



noiseless 

low noise 

high noise 

RMS 

on a 

0.2% 

38.6% 

127.3% 

RMS 

on / 

0.13% 

18.7% 

55.1% 


Table 1. Relative Mean-Square errors on a and / at 8 iterations corresponding to the 
plots displayed on Fig. 6. 


(ii) As may be seen on Fig. 6, althongh strong additive noise impacts the 
reconstrnctions badly, one may notice on the cnt plots that the oscillations on 
reconstrnctions average abont each constant valne, leading us to believe that a 
penalization term (e.g. total variation norm) favoring piecewise constant functions 
would re-establish good convergence. Additionally, noise in data may make the 
algorithm give negative values to both a and /, although both quantities are 
physically nonnegative. This may be avoided by introducing at each iteration a 
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0 0.5 1 0 0.5 1 -1 0 1 2 



0 0.5 1 0 0.5 1 - 0.5 0 0.5 1 1.5 


Figure 6. Reconstructed a (top row) and / (bottom row) after convergence. Left to 
right: noiseless, low noise, high noise, cut plots at {a; = 0}. 

projection step onto nonnegative functions (i.e. of the form a{x) = max(a(a;), 0)), 
at the cost of losing the first property of averaging around the correct constant 
values. Finding an algorithm taking additive noise into account while respecting 
physics-based criteria appropriately will be the object of future work. 
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